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OfT-equilibrium dynamics of a three-dimensional lattice model with nearest- and next nearest- 
neighbors exclusions is studied. At equilibrium, the model undergoes a first-order fluid-solid tran- 
sition. Non-equilibrium filling, through random sequential adsorption with diffusion, creates amor- 
phous structures and terminates at a disordered state with random closest packing density that lies 
in the equilibrium solid regime. The approach towards random closest packing is characterized by 
two distinct power-law regimes, reflecting the formation of small densely packed grains in the long 
time regime of the filling process. We then study the fixed-density relaxation of these amorphous 
structures towards the solid phase. The route to crystallization is shown to deviate from the simple 
grain growth proposed by classical nucleation theory. Our measurements suggest that relaxation is 
driven mainly by coalescence of neighboring crystallized grains which exist in the initial amorphous 
state. 



INTRODUCTION 

Random sequential adsorption (RSA) processes are 
well-studied examples of off-equilibrium dynamics [T] . In 
these processes, particles are sequentially introduced into 
a system at random locations, and get adsorbed if this 
location is compatible with the other particles. For mod- 
els with inter-particle hard-core exclusion, the process 
ends when all deposition options are exhausted. Typi- 
cally, RSA processes terminate at a limiting density far 
lower than pcp, the closest packing density of the re- 
spective system. For example, random deposition of 
dimers on a one-dimensional (ID) lattice creates sites 
whose two nearest neighbors are filled. These sites can 
never be filled and the process jams with terminal density 

e., ^ 0.864664 ig. 

When adsorbed particles are allowed to diffuse, 
jammed states can be relaxed by particle movement 
(RSAD models). Jamming density in RSAD models is 
naturally higher and in some cases closest packing den- 
sity is reached. In some cases, such as the ID lattice 
RSAD and in the square lattice exclusion models Nl 
and N2 (particles deposited on a square lattice, where 
each particle excludes deposition on its four nearest- 
neighbors (Nl) or up to second nearest neighbors (N2)), 
the density converges to pcp via a power law 

[3H5]. In contrast, in the 2D N3 lattice model (exclu- 
sion up to third nearest neighbors) the density converges 
via & power law to a limiting density lower than 
Pcp J which is termed the random closest packing den- 
sity Prop = 0.171626 < pcp = 0.2 0. 

Random closest packing states are characterized by the 
absence of mobile particles (but a few rattlers). Thus, 
diffusive relaxation from these off-equilibrium states is 
impossible. However, off-equilibrium diffusive relaxation 
can be studied by terminating RSAD filling process at 
a predetermined density, lower than the random closest 
packing one, and following relaxation dynamics towards 



crystallization. Better characterization of the crystalliza- 
tion process is of interest (see e.g. [3-^), as some of its 
aspects, most notably the crystallization rate, are not 
well explained by standard nucleation theory In ad- 
dition, recent works have highlighted the possibility that 
while classical nucleation theory describes well the nu- 
cleation at concentrations slightly above the coexistence 
region, highly saturated super-cooled fluids undergo dif- 
ferent mechanisms of nucleation. These include spinodal 
nucleation and a third mechanism, yet to be fully under- 
stood, which is characterized by small displacements of 
the particles PUIITT] . 

Here we study a three-dimensional (3D) model on 
the cubic lattice, where the inter-particle interaction is 
limited to hard-core exclusion up to the second-nearest 
neighbors, henceforth termed the 3DN2 model. This can 
be thought of as a system of hard particles of volume 
4 sites, whose closest packing configuration is a body- 
centered cubic (BCC) arrangement with a lattice con- 
stant equal to 2 lattice sites (Fig. [T]) . The closest packing 
particle density is therefore pcp = 0.25. In the following 
we measure all densities relative to this maximal parti- 
cle density (n = p/pcp), be. we use volume density. 
We start by exploring the equilibrium behavior of the 
3DN2 model, and show that it undergoes a first-order 
phase transition, with a density gap between Uf = 0.413 
and Us — 0.519. We then look at the RSAD dynam- 
ics, and show that unlike the 2D N3 model, where the 
random closest packing density lies within the density 
gap, the RSAD dynamics of 3DN2 terminates at den- 
sity njicp = 0.81, well within the solid phase. Thus, 
RSAD dynamics enables us to compress the system deep 
into the super-cooled regime, while keeping the system 
microscopically amorphous. 

Finally, we study the diffusive relaxation of the amor- 
phous state as a function of density. The study of the 
off-equilibrium relaxation of disordered super-cooled con- 
figurations towards the equilibrium solid has attracted 
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much interest, mainly in the context of super-cooled 
glass- formers near the glass transition [l2l - fT5] . In con- 
trast, here this relaxation process leads into the equilib- 
rium crystal. Yet, it turns out that several characteristics 
of the dynamics are reminiscent of glassy relaxation. In 
fact, it has been argued that in a number of systems 
an ideal glass transition could be related to the kinetic 
spinodal [101 fT6H20] , in the sense that spontaneous crys- 
tallization of many stable micro-crystallites could lead to 
a kinetic arrest due to the competition between the dif- 
ferent crystalline orders. In comparison, the 3DN2 sys- 
tem presents no such competition, and thus the kinetic 
spinodal is manifested only in the form of a modified so- 
lidification process. Indeed, we find that the route to 
crystallization in the high densities, well above the equi- 
librium transition density, is characterized by most parti- 
cles moving as little as needed to shift to the sublattice on 
which the system crystallizes, in agreement with recent 
works on glassy hard spheres system [llj . Furthermore, 
we observe that a leading mechanism is the alignment of 
neighboring micro-crystallites, which exist in the amor- 
phous initial state, against each other. 
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FIG. 1. A hard core particle which enforces the model's ex- 
clusion rules (left). The central cube is one cubic-lattice site 
and the whole particle spans 2 lattice sites in each direction. 
Several configurations are shown in order to demonstrate the 
allowed nearest neighbors and the closest packing configura- 
tion - a BCC lattice (right). 



THERMODYNAMICS 

To characterize the model's thermodynamics, we used 
Monte-Carlo (MC) simulations of a grand canonical en- 
semble. In each MC step, particles are allowed to be 
adsorbed on the lattice (while adhering to the exclusion 
rules) or desorbed from it. The probability of adsorb- 
ing/desorbing a particle is governed by the chemical po- 
tential fi. We started the simulations with an empty 
lattice and raised the activity z = e^'^ monotonically 
(compressing, or cooling), equilibrating the system af- 
ter each raise. At the critical regime, the distribution 
of the density for a finite system (at a constant z) is bi- 
modal (Figure [2|^ a)), reflecting the transitions between 
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FIG. 2. (Color online) The phase transition, (a) A histogram 
of the density at constant z ~ Zc showing the system divides 
its time between two, high (solid) density and low (liquid) 
density, states, meaning the transition is of first order, (b) 
The average density in the vicinity of the transition as mea- 
sured in lattices of different linear size. The inset shows the 
finite size scaling of Zc according to the locations of the peaks 
in the compressibility (dashed line) and of the dips in Binder's 
cumulant (Eq. [2| (solid line), (c) and (d) are the measured 
K.{z) and Ul{z), respectively, and (in the insets) the finite size 
scaling of their critical values. 
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two phases of different density - higher density sohd and 
lower density fluid. Thus, the transition is a first order 
one. 

We used finite size scaling, with lattices of sizes L = 
12,14,16... 30, to extrapolate the transition parameters 
at the thermodynamic limit [H]. This analysis (Figure 
[2]) verified that the transition is indeed first order. The 
compressibility, 



(1) 



exhibits a rounded peak at the transition, whose height is 
proportional to the volume, as expected for a first-order 
transition. Binder's cumulant, the forth order moment 
of the form 
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also conforms to the first-order behavior - it is constant 
far from the transition, and near it it has a dip whose 
height is proportional to [53]. The location of the 
transition, can be extracted from either the positions of 
k's maxima or of J7l's minima. Both converge (the latter 
faster than the first) towards a value of Zc = 1.7010(5). 
The limiting densities of the liquid and solid phases are 
those around which the distribution of the density peaks, 
Uf = 0.4131 ± 0.003 and = 0.5195 ± 0.0005 accord- 
ingly, in agreement with a previous thermodynamic study 
|22j . Here too, the finite size correction exhibits the non- 
anomalous 0{L~^) behavior. 

APPROACH TO RANDOM CLOSEST PACKING 

We have studied the RSAD dynamics |5J of the 3DN2 
system in extensive MC simulations, using lattices of lin- 
ear size between 64 and 256 (over 3 million particles). 
RSAD dynamics is equivalent to a grand canonical dy- 
namics with an infinite chemical potential, resulting in 
particle deposition whenever possible, and no desorption. 
As the only energy scale in the problem is the chemical 
potential, infinite chemical potential is equivalent to an 
infinitely fast cooling. First, we performed random se- 
quential deposition on the (initially empty) lattice. This 
filling process happens at time zero and terminates at 
density 0.49362 ± 0.00004, in the equilibrium coexistence 
regime. When all adsorption possibilities are exhausted, 
the system is let to diffuse, i.e. particles attempt to move 
to a randomly selected nearest neighbor site. Whenever 
diffusion creates a space large enough for a particle to be 
added without violating the exclusion rules, a particle is 
immediately placed there. A rejection-free algorithm was 
used to reduce run-time and allow simulation of larger 
lattices. 

Measurements of volume density (Fig. [3]) clearly show 
that the process jams at a limiting density rircp — 
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FIG. 3. Random sequential adsorption with diffusion 
(RSAD). The volume density is plotted vs. time on a loga- 
rithmic scale. The results are averaged over 1500 realizations, 
for lattice of size L = 256 (finite size effects were minute, thus 
we present only the results for L = 256). Relative standard 
error is less than 0.0005 for all t. 
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FIG. 4. (Color online) Density vs. time along the RSAD 
process. The residual density (distance, in density, to the ter- 
mination random closest packing density) is plotted, in order 
to highlight the two power-law regimes: the earlier with a 
residual density proportional to t~^^^ (dashed line), and the 
later one exhibiting a behavior (dotted line). 



0.81095 ± 0.00004, lower than the closest packing value, 
but much higher than the transition region. The long- 
time approach towards rircp follows the asymptotics 
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-A/t (see Fig. HI). However, at shorter times. 



one observes a regime in which the residual vacancies den- 
sity follows a t~^^^ rule. The long-time RSAD dynamics 
can be described by a model of holes diffusion (facili- 
tated by diffusion of the neighboring particles) |24j. The 
vacancies in the lattice can be relaxed by particle depo- 
sition only when multiple holes colocalize, which can be 
described as an annihilation process from the holes' per- 
spective. The simplified diffusion-annihilation descrip- 
tion predicts that the vacancies density should follow a 
t^'^/^ law in 3D {t~^ in 2D), according to the probability 
of two diffusing holes to meet (or the return probability 
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of a random walk). Here we find that the intermediate- 
time behavior agrees with hole diffusion in 3D while the 
long-time 1/t behavior suggests that hole diffusion is ef- 
fectively 2D. Analyzing the microscopic processes, we 
found that this effective dimensional reduction results 
from the formation of stable, perfectly packed, clusters. 
The holes left in the systems are trapped between clus- 
ters and perform 2D diffusive motion along the bound- 
aries of the clusters, resulting in a slower density gain, 
rircp — n{t) ^ t~^, associated with the effectively 2D dif- 
fusion. A similar phenomenon was observed for the 2D 
Nl and N2 models, where long-time RSAD relaxation 
was governed by ID processes f5]. 

The existence of packed clusters at long times is cor- 
roborated by measurements of the static density-density 
correlation function 

(p(r)p(r + R))-(p(r))^ 
^^^-^tW = ^(^) ' (3) 

where 

is the total (time independent) number density. Taken 
along the principal axes and main diagonals, the cor- 
relation function is characterized by the peaks and dips 
manifesting the closest packing BCC order in short range: 
peaks for even integers and dips for odd integers (princi- 
pal axes) as well as peaks for multiples of -s/S (diagonals). 
The closely-packed clusters grow in the short time regime 
until they are stopped by neighboring clusters residing in 
other sublattices. The decay of the correlation function 
shows that the typical size of the clusters towards the end 
of the filling process is about 5 sites. The clusters reach 
their maximum size when the density is a little over 0.8, 
coinciding with the later time regime when density gain 
starts converging through the rule (Fig. [s]). 

RELAXATION DYNAMICS 

As the random closest packing density lies deep within 
the solid phase, RSAD dynamics may be used to cre- 
ate a disordered configurations at densities deep within 
the equilibrium solid regime. In order to study off- 
equilibrium relaxation of the super-cooled fluid, we there- 
fore employ RSAD dynamics to fill the system up to 
some fixed density, and then let the system relax dif- 
fusionally, keeping the number of particles fixed. An 
important characterization of the dynamics is given by 
particles' mobility. Off-equilibrium systems are often 
characterized by a sub-diffusive behavior. For example, 
some glassy systems exhibit caged dynamics where parti- 
cles are trapped in small regions, within which diffusion 
takes place. Such systems display a short-time diffusion 




FIG. 5. (Color online) The static correlation function gstat(R.) 
at different times along the filling process, taken along the 
principal axes and main diagonals and plotted against i?=|R| . 
The positive peaks (at R's fitting the closest packing config- 
uration) decay exponentially with distance. The correlation 
length reaches A ~ 5.3 sites at pRCP (regardless of lattice 
size), as shown by the fit to a pure exponent (dashed line). 

within the cages, followed by a period characterized by 
sub-diffusion as a result of the localization of the parti- 
cles. The transition rates between cages determines the 
long-range, long-time, mobility and could be either dif- 
fusive or sub-diffusive depending on the nature of the 
barriers between the cages. Here, the average squared 
displacement displays a similar behavior, as shown in 
Fig. [6] only that the later stages are inter-mixed with 
the solidification process itself. The system starts at a 
disordered state so that the particles are caged and lo- 
calized - free brownian movement is only measured at 
very short times, of order of one time unit, meaning the 
average cage size is of order a single lattice site. Then, af- 
ter some (density-dependent) time, structural relaxation 
begins with movement between cages, and the diffusion 
measured is an anomalous diffusion resulting from the 
distribution of cages' size. The higher the density, the 
lower is the exponent characterizing the sub-diffusive mo- 
tion (see Fig. |6|. Finally, as the system crystallizes, par- 
ticles get localized permanently and mobility diminishes. 

The evolution of the relaxation process can be quanti- 
tatively studies through the density-density autocorrela- 
tion function 



C{t) = -^^—^ \^-Y,p{n,0)p{r,,t)-p'j , (5) 

where runs over all V sites, and p{vi,t) is one (zero) 
if at time t there is (not) a particle at site r,. This cor- 
relation function (or rather its long-time limit) is often 
used as the Edwards- Anderson order parameter in glassy 
materials. Beyond the glass transition, finite correlations 
are expected even for infinitely long times. Apparently, 
the model exhibits the two-step relaxation that is typical 
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FIG. 6. (Color online) The mean squared displacement aver- 
aged over all particles in ten realizations of a system of size 
L — 128, as a function of density. The log-log plot high- 
lights the diffusive regimes: in times of order of a single time 
step, the particles perform free diffusion within cages. As the 
cages' sizes are of order of 1 site, this is followed by a sub diffu- 
sive stage, where the particles are still mostly bound to their 
cages, and vibrate around their origins. At later times, the 
diffusivity is dominated by the movement of particles between 
cages and greatly depends on density. For 10"^ < t < 10^, we 
see that at densities within the transition gap the behavior is 
again diffusive (with a diffusion constant 100 times smaller 
than that of the first regime) , but systems that were super- 
cooled into the solid regime exhibit sub-diffusive behavior, 
characterized by a power- law < >= Dt" with a < 1. Fi- 
nally, for long times, diffusion is diminished as the systems 
goes through solidification (The longer times are inaccessible 
in the gap densities as some of the particles already reach the 
limit of the finite system. 



of glassy systems (Fig. [?]). Due to the aforementioned 
caging phenomenon, there is very little relaxation in the 
short time regime as the particles mostly move back and 
forth within their cages. This leads to the formation of 
a plateau in the autocorrelation function, known as the 
/3-regime. As density increases, the plateau extends to 
longer times. It is followed by the structural relaxation 
in the a-regime, where cages are broken and restructured, 
and the correlation function decays to zero. While some 
of the simulations ended having finite C{t) a detailed 
finite size scaling shows that this is a finite size effect. 
The a-decay takes the form of a stretched exponential, 
expressing the spatially heterogeneous character of the 
system. 

Despite the apparent similarities to the behavior of 
glassy systems, it should be stressed that the long-time 
relaxation observed here is qualitatively different from 
the alpha-phase relaxation seen in glassy materials. As 
indicated in Fig. [7j the long time relaxation overlaps 
with the solidification of the system. The nature of this 
process and the reason it is accompanied by spatially 
long-ranged relaxation can be understood looking at the 
microscopic picture. 



FIG. 7. (Color online) C{t) in an L = 128 Lattice, for den- 
sities ranging from 0.45 to 0.8. The results are an average 
over several hundred realizations (there are less realizations 
in the long time regime, hence the larger error bars). Black 
line indicating the time at which half of the particles are on 
the same sublattice. 
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FIG. 8. (Color online) Particles' displacements in a high- 
density (n = 0.75) system at a time it has already crystallized 
(almost all particles are lying on the same BCC sublattice). 
Most particles have moved less than \/3 sites. The inset shows 
the displacement correlation length as estimated by the decay 
of the displacement correlation function to 0.2. Displacement 
correlation is short ranged even for times longer than crystal- 
lization time. Data presented for a system of size L — 128 
and density n = 0.75 (over 100 realizations), and is typical 
for other high densities and reasonable correlation cutoffs. 



MICROSCOPIC PICTURE OF THE 
CRYSTALLIZATION DYNAMICS 

Detailed analysis of microscopic processes underlying 
the dynamics reveals some interesting features of the re- 
laxation process. We focus on the total displacements of 
particles along the process (Fig. [s]), and find that in the 
high-density regime a considerable number of the par- 
ticles move only a distance of 1 to -s/S sites away from 
their t — location, i.e. up to one shift along each of 
the axes. Moreover, the particles' movement distance is 
highly correlated between neighboring particles, i.e. dif- 
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ferent regions in the model have their own characteristic 
distance. The sizes of these clusters can be quantified by 
the spatial correlations of the displacement 

Cc(r,t) = (i?(0,t)i?(r,t))-(i?)2, (6) 

which exhibit short-range correlations. In order to ob- 
tain an estimate for the length scale characterizing these 
displacement correlation we look at the distance at which 
CD(r, t) decays to an arbitraty cutoff value (0.2). The in- 
set of FigjS] presents the time-dependence of this length 
scale. Not only are the estimated correlation lengths 
short-ranged compared to the static correlation measured 
in t=0, but they also do not grow over time. It is im- 
portant to notice that even though the exact value of the 
length scale is as arbitrary as the threshold used, these 
two observations are not. 

The above observations are inconsistent with simple 
nucleation theory picture of crystallization, according to 
which growing nucleation seeds are continuously fed by 
their fluid environment. Instead, we observe the align- 
ment of micro-crystallites (present at t = 0) to be the 
main mechanism leading to crystallization. The par- 
ticles forming the boundaries of the clusters are effec- 
tively caged by the surrounding clusters, so that their 
movement is somewhat restricted (resulting in the sub- 
diffusive behavior described above). This microscopic 
picture fits with the global dynamic measurements de- 
scribed above: the displacement correlation length ob- 
served above reflect nothing but the size of the clusters 
moving cooperatively to re-align themselves against each 
other. The crystallization mechanism described here, 
characterized by minute displacements of particles along 
the process, is reminiscent of a similar phenomenon re- 
cently observed for the hard spheres glassy system [TT] . 
It is believed that the transition to this regime has to do 
with the density being so high (or the super-cooled fluid 
being over-saturated) such that the kinetic spinodal is 
reached. Under these conditions, micro-crystallites of 
any size are stable. Thus, virtually all of the micro- 
crystallites naturally occurring in the initial amorphous 
state survive the relaxation process, and crystallization 
is achieved by re-alignment of these against each other. 
This highlights another possible connection to the glassy 
behavior. It was suggested that in some systems an ideal 
glass transition occurs at the kinetic spinodal, as micro- 
crystallites having competing crystal order are all stable 
and thus no global order may be achieved [l6U20j . The 
fate of a deeply super-cooled fluid is then determined by 
the nature of the micro-crystallites occurring when ki- 
netic spinodal is hit. If they are easily alignable, one 
would end up having a solid, while when they are not 
easily alignable, a stable amorphous glass would emerge. 

Finally, we comment about the crystallization time. 
We define the crystallization time as the time (measured 
from the end of the RSAD process) till half of the par- 
ticles occupy the same sublattice (The equilibrium oc- 
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FIG. 9. (Color online) The time of solidification, defined as 
the first time at which more than half of the particles are 
on the same sublattice, averaged over 77 realizations, as a 
function of system size and density. Realizations in which two 
or more macroscopic domains are formed are averaged only if 
one of the domains contains more than half of all particles. 



cupation fraction of the leading sublattice depends only 
weakly on the density, changing from 0.9 at the onset of 
crystallinity to 0.96 for n = 0.6 and approaching unity 
for high densities). Measuring this time as a function of 
density, one observe a non- monotonous picture (Fig. |9|. 
As density increases, the preference of the solid phase 
over the disordered phases increases. At the same time, 
the dynamics becomes increasingly sluggish. The result 
of these two effects combined is the non-monotonous be- 
havior. The long relaxation time in high-density systems 
also stems from the fact that they are prone to the forma- 
tion of two or more macroscopic domains which require 
a long time to align with one another. The increase in 
the average solidification time is enhanced by the finite 
size effects - the larger the lattice, the longer it takes for 
the system to solidify. 



CONCLUSION 

The 3DN2 model presented here exhibits a non triv- 
ial RCP dynamics. The density gain in the long-time 
regime, controlled by holes' diffusion, is going through 
two distinguished regimes corresponding to an unordered 
dense liquid state, followed by a state of many micro- 
crystallites in which the particles mainly perform 2D dif- 
fusion on the clusters' boundaries. Keeping the density 
fixed, the system relaxes to a solid. In the high-density 
regime, beyond the kinetic spinodal, even small micro- 
crystallites are stable. Relaxation is then brought about 
through alignment of the small solid clusters against each 
other. It then follows that the dynamics of the system 
is heterogeneous, similarly to super-cooled glass- formers. 
The crystallization process, is dramatically different than 
the one described by classical nucleation theory. 
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